/// [Optimum polynomial](https://projecteuler.net/problem=101)
fn main() {
    for x in &[
        vec![0.0, 0.0, 0.0, 1.0],
        vec![1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0],
    ] {
        let base = mat_org(x.len());
        if let Ok(data) = mat_mul(&base, &x) {
            let mut sum = 0;
            for n in 1..x.len() {
                if let Some(src) = data.get(0..n) {
                    if let Ok(rst) = optimum_polynomial(src) {
                        sum += rst.round() as i64;
                    }
                }
            }
            println!("{}th degree polynomial generating function:\t{}", x.len(), sum);
            let mut sum = 0;
            for n in 1..x.len() {
                if let Some(src) = data.get(0..n) {
                    sum += lagrange_interpolation_polynomial(src).round() as i64;
                }
            }
            println!("{}th degree polynomial generating function:\t{}", x.len(), sum);
        }
    }
}
/// 拉格朗日插值法
fn lagrange_interpolation_polynomial(src: &[f64]) -> f64 {
    let len = src.len() as i32;
    let mut sum = 0.0;
    for i in 0..len {
        let mut mul = 1.0;
        for j in 0..len {
            if i != j {
                mul *= (len - j) as f64 / (i - j) as f64;
            }
        }
        sum += mul * src[i as usize];
    }
    sum
}
/// 矩阵运算，逆矩阵
fn optimum_polynomial(src: &[f64]) -> Result<f64, &'static str> {
    match mat_inv(&mat_org(src.len())) {
        Ok(base) => {
            let x: Vec<f64> = src.to_vec();
            match mat_mul(&base, &x) {
                Ok(rst) => {
                    let mut sum = 0.0;
                    for i in 0..rst.len() {
                        sum += rst[i].round() * ((rst.len() as i64 + 1).pow(i as u32)) as f64;
                    }
                    Ok(sum)
                }
                Err(e) => Err(e),
            }
        }
        Err(e) => Err(e),
    }
}
/// 获取 x.pow(n) 矩阵
fn mat_org(n: usize) -> Vec<Vec<f64>> {
    let mut mat: Vec<Vec<f64>> = vec![];
    for i in 1..=n {
        let mut tmp: Vec<f64> = vec![];
        for j in 0..n {
            tmp.push(i.pow(j as u32) as f64);
        }
        mat.push(tmp);
    }
    mat
}
/// 求解逆矩阵
fn mat_inv(src: &Vec<Vec<f64>>) -> Result<Vec<Vec<f64>>, &'static str> {
    if src.len() > 0 && src.len() == src[0].len() {
        let mut inv: Vec<Vec<f64>> = vec![];
        for i in 0..src.len() {
            let mut tmp: Vec<f64> = vec![0.0; src.len()];
            tmp[i] = 1.;
            inv.push(tmp);
        }
        let mut src = src.to_vec();
        let mut stat = true;
        'outer: for i in 0..(src.len() - 1) {
            let j = i + 1;
            let a = src[j][i] / src[i][i];
            let b = src[j][j] - src[i][j] * a;
            for k in 0..src.len() {
                if b == 0. {
                    stat = false;
                    break 'outer;
                }
                src[j][k] = (src[j][k] - src[i][k] * a) / b;
                inv[j][k] = (inv[j][k] - inv[i][k] * a) / b;
            }
            for j in (i + 2)..src.len() {
                let a = src[j][i] / src[i][i];
                for k in 0..src.len() {
                    src[j][k] = src[j][k] - src[i][k] * a;
                    inv[j][k] = inv[j][k] - inv[i][k] * a;
                }
            }
        }
        for i in (1..src.len()).rev() {
            for j in (0..i).rev() {
                let a = src[j][i];
                src[j][i] = src[j][i] - src[i][i] * a;
                for k in 0..src.len() {
                    inv[j][k] = inv[j][k] - inv[i][k] * a;
                }
            }
        }
        if stat {
            Ok(inv)
        } else {
            Err("no solution")
        }
    } else {
        Err("src.rows != src.cols")
    }
}
/// 矩阵乘法
fn mat_mul(two: &Vec<Vec<f64>>, one: &Vec<f64>) -> Result<Vec<f64>, &'static str> {
    if two.len() > 0 && two[0].len() > 0 && one.len() == two[0].len() {
        let mut dst: Vec<f64> = vec![];
        for rows in two {
            let mut sum = 0.0;
            for i in 0..one.len() {
                sum += rows[i] * one[i];
            }
            dst.push(sum);
        }
        Ok(dst)
    } else {
        Err("base.rows != x.cols")
    }
}
